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Abstract 

This paper advocates use of the atomic orbitals which have direct physical interpretation, 
i.e. Coulomb Sturmians and hydrogen-like orbitals. They are exponential type orbitals (ETOs). 
Their radial nodes are shown to be essential in obtaining accurate nuclear shielding tensors for 
NMR work. The present work builds on a 2003 French PhD and many numerical results were 
published by 2007. The improvements in this paper are noteworthy, the key being the actual 
basis function choice. 

Until 2008, their products on different atoms were difficult to manipulate for the evalua- 
tion of two-electron integrals. Coulomb resolutions provide an excellent approximation that 
reduces these integrals to a sum of one-electron overlap-like integral products that each involve 
orbitals on at most two centers. Such two-center integrals are separable in prolate spheroidal 
co-ordinates. They are thus readily evaluated. Only these integrals need to be re-evaluated to 
change basis functions. 

In this paper, a review of the translation procedures for Slater type orbitals (STO) and 
for Coulomb Sturmians follows that of the more recent application to ETOs of a particularly 
convenient Coulomb resolution. 



Keywords: Coulomb Sturmian basis, nodal structure, Coulomb resolutions, ab initio quan- 
tum chemistry. 
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1 Introduction 



The criteria for choice between gaussian and exponential basis sets for molecules do not seem obvious 
at present. In fact, it appears to be constructive to regard them as being complementary, depending 
on the specific physical property required from molecular electronic structure calculations. 
The present work describes a breakthrough in two-electron integral calculations, as a result of 
Coulomb operator resolutions. This is particularly significant in that it eliminates the arduous 
orbital translations which were necessary until now for exponential type orbitals. The bottleneck 
has been eliminated from evaluation of three- and four- center integrals over Slater type orbitals 
and related basis functions. 

The two-center integrals are replaced by sums of overlap-like one-electron integrals. This implies a 
speed-up for all basis sets, including gaussians. The improvement is most spectacular for exponen- 
tial type orbitals. A change of basis set is also facilitated as only these one-electron integrals need 
to be changed. The gaussian and exponential type orbital basis sets are, therefore interchangeable 
in a given program. The timings of exponential type orbital calculations are no longer signifi- 
cantly greater than for a gaussian basis, when a given accuracy is sought for molecular electronic 
properties. 

Atomic orbitals are physically meaningful one-electron atom eigenfunctions for the Schrodinger 
equation. This gives well-known analytical expressions: hydrogen-like orbitals. 
Boundary conditions allow the principal quantum number n to be identified as the order of the 
polynomial factor in the radial variable. It must therefore be positive and finite. It is also defined 
such that n — I — I is greater than or equal to 0. This gives the number of zeros of the polynomial 
(radial nodes). Here, I = or a positive integer, which defines the angular factor of the orbital, (i.e. 
a spherical harmonic, or, more rarely, its Cartesian equivalent) The number n gives the energy of 
the one-electron atomic bound states. Frequently, basis set studies focus on the radial factor. That 
is, for our present purposes, the angular factor can be assumed sufficiently defined as a spherical 
harmonic. 

The key issue is whether to choose basis sets with exponential or gaussian asymptotic factors. 
Certain physical properties, such as NMR shielding tensor calculations directly involve the nu- 
clear cusp and correct treatment of radial nodes, which indicates that basis sets such as Coulomb 
Sturmians are better suited to their evaluation than gaussians [HSJ GUJ SB] • 

There is also evidence to suggest that CI expansions converge in smaller exponential basis sets 
compared to gaussians \A5\ WQ . Benchmark overlap similarity work is available |46t 

2 Wave-function quality 

The following quantity: 



is used to test wave-function quality. It is smooth, to varying degrees, in different basis sets. 
Atomic positions must give cusps. The importance for stable and accurate kinetic energy terms, 
particularly in DFT. 

Much molecular quantum chemistry is carried out using gaussian basis sets and they are indeed 
convenient and lead to rapid calculations. The essential advantage they had over exponential basis 
sets was the simple product theorem for gaussians on two different atomic centers. This allows all 
the two-electron integrals, including three- and four- center terms to be expressed as single-center 
two-electron integrals. 
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The corresponding relationship for exponential type orbitals generally led to infinite sums and the 
time required, particularly for four-center integrals could often become prohibitive. 
Recent work by Gill has, nevertheless been used to speed up all three and four center integral 
evaluation, regardless of basis using the resolution of the Coulomb operator [15\ \W[ I36| . This 
work by Gill is used here to reduce the three- and four- center two-electron integrals to a sum of 
products of overlap-like (one-electron) integrals, basically two-centered. This algorithm was coded 
in a Slater type orbital (STO) basis within the framework of the STOP package [I] (in fortran) 
during summer 2008. Note, however, that other exponential or gaussian basis sets can readily be 
used. The set of one-electron overlap-like auxiliary integrals is the only calculation that needs to 
be re-done. They may be re-evaluated for the basis set that the user selects for a given application. 
This procedure makes the approach highly versatile, since a change of basis set requires relatively 
few simple new evaluations. A modular or object oriented program is being designed to do this 
efficiently [361 SSI H] ■ 

The present article gives illustrative test results on molecular systems e.g. the H2 dimer. 

The layout is as follows: the review begins with a brief recap of basis sets and programming 
strategy in the next two sections. Atom pairs are the physical entity used for integral evaluation, 
both in the Poisson equation technique and the Coulomb resolution. Two sections are devoted 
to these progressively more powerful techniques which both reduce two-electron to one-electron 
integrals. The overlaps required for the Coulomb resolution differ by a potential factor from orbital 
overlaps. Their evaluation is nevertheless analytic, using well-known techniques summarized in the 
subsequent section. Finally, to illustrate what can be gained by eliminating orbital translations, 
the translation of Slater type orbitals is reviewed briefly, from recent work on BCLFs. A few 
numerical results are given on the dimer of molecular hydrogen which show progressive speed-up 
particularly for the Coulomb resolution given a pre-selected accuracy, which proves sufficient to 
provide satisfactory confirmation of experimental work on this dimer. 

3 Basis sets 

Although the majority of electronic quantum chemistry uses gaussian expansions of atomic or- 
bitals |191l20j. the present work uses exponential type orbital (ETO) basis sets which satisfy Kato's 
conditions for atomic orbitals: they possess a cusp at the nucleus and decay exponentially at 
long distances from it [58l \59\ [60]. It updates a tradition beginning around 1970 and detailed 
elsewhere [gSlEIHEQlEZllZIllZflESESEElEnE^ 

Two types of ETO are considered here: Slater type orbitals (STOs) [56tl57| and Coulomb Sturmians 
and their generalisation, which may be written as a finite combination thereof |102| . Otherwise, 
STOs may be treated as multiple zeta basis functions in a similar way to the approach used with 
gaussian functions. 

Many exponential type functions exist [102] . Preferential use of Sturmian and related functions 
with similar radial nodes is discussed [36j . 

Coulomb Sturmians have the advantage of constituting a complete set without continuum states 
because they are eigenfunctions of a Sturm-Liouville equation involving the nuclear attraction 
potential i.e. the differential equation below. 

v^;mo = [/3 2 -^]^(/?,r). 

The exponential factor of Coulomb Sturmians; e~^ T has an arbitrary screening parameter j3. In 
the special case when /3 = C/n with n the principal quantum number and £ the Slater exponent, 
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we obtain hydrogen-like functions, which do not span the same space and require inclusion of 
continuum states to form a complete set |102| . Hydrogen-like functions are, however well known as 
atomic orbitals: the radial factor contains the associated Laguerre polynomial of order 21 + 1 with 
suffix n — l — 1 and the exponential e~'* r / n as indicated above. The angular factor is just a spherical 
harmonic of order I. These functions are ortho- normal. The optimal values of the f3 parameters 
may be determined analytically by setting up secular equations which make use of the fact that 
the Sturmian eigenfunctions also orthogonalise the nuclear attraction potential, as developed by 
Avery [71]. 

J S£(r,0,0)3&M,0)y = <W«W- 

These ortho-normal functions are further generalised by varying a from the Coulomb Sturmian value 
of 1. In such a case, the basis remains ortho- normal and othogonalises a/r a . This eliminates the 
r 2 term, arising for quadrupole moments when a = —2, thus confirming the very recent numerical 
observations in the Guseinov group |107| . Similarly, it would be expected that a = — 1 ETOs 
would constitute the optimal basis for magnetic dipole integrals of NMR shielding. Furthermore, 
a negative value of a will not modify the number of radial nodes: the functions will simply breath. 
Definitions: the generalised exponential functions constitute finite complete orthonormal sets. 
Their expression is as follows: 



Xnlm (r) = [^] N nl L^_~ a _ a (2 Cr) r'e^Yf <t>) (3-1) 
y In 

Here, N is the normalisation constant previously obtained for Coulomb Sturmians, L is the associ- 
ated Laguerre polynomial of order 21 + 2 — a with suffix n — I — 1 — a (recall that a = 1 defines the 
Coulomb Sturmians. 

Define a variable including the screening constant: 

x = 2 Cr 

Subsequently, rewriting the norm as N(a) n [ and introducing p = 21 + 2 — a and q = n + l + l — a 
gives the simplified expression for the generalised orthonormal basis sets of ETO, used by Guseinov. 



XnlmW = N(a) nl L? (x) r l e- x / 2 Yr(8, <f>) (3.2) 

In past applications, no obvious advantage has been evidenced for the functions with negative a 
indices over the well-known Slater type orbitals, Coulomb Sturmians (a = 1) or Shull-Loewdin 
functions (a = 0). In fact, the infinite series arising when Hartree-Fock two-electron integrals that 
do not possess closed forms (three and four center terms) are evaluated converge much more slowly 
when the negative a functions are used. This has recently also proven to be the case of a set of 
electric field integrals [108] . 

This paper records the precedent of electric quadrupole integrals, already published by Guseinov's 
colleagues, where the negative a basis converges as well as (if not better than) the STO [107] and 
presents a new application to the dipole integrals in the NMR experimental setup. 
The investigations are extended to comparisons with previous work on the nuclear dipole integrals 
that are so important to the evaluation of nuclear shielding tensors and NMR chemical shifts. 
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Furthermore, these nuclear magnetic dipole integrals are closely related to the one-electron nuclear 
attraction integral, required in all Hartree-Fock and DFT work. 

Alternative ETOs would be Slater type orbitals and B-functions with their simple Fourier trans- 
forms. Strictly, they should be combined as linear combinations to form hydrogen-like or, better, 
Sturmian basis sets prior to use. 

STOs allow us to use routines from the STOP package [55] directly, whereas Coulomb Sturmians 
still require some coding. The relationship to STOs is used to carry out calculations over a Coulomb 
Sturmian basis with STOP until the complete Sturmian code is available. The present state-of-the- 
art algorithms require at most twice as long long per integral than GTO codes but the CI converges 
with fewer functions and the integrals may be evaluated after gaussian expansion or expressed as 
overlaps to obtain speed up |63j. 



After a suitably accurate electron density has been obtained for the optimized geometry over a 
Coulomb Sturmian basis set, the second order perturbation defining the nuclear shielding tensor 
should be evaluated in a Coupled perturbed Hartree Fock scheme. 

The integrals involved may conveniently be evaluated using B-functions with linear combinations 
giving the Coulomb Sturmians. 

S nl (r)-(2a) 2l + 1)ll t\(l + S/2) t Bt+l ^ ] 

The techniques exploit properties of Fourier transforms of the integrand. 

Note that either HF or DFT can serve as zero order for the present nuclear shielding tensor calcu- 
lation over ETOs. 

A full ab initio B-function code including nuclear shielding tensor work is expected to be complete 
shortly. 

Some tests show that Slater type orbitals (STO) or B-functions (BTO) are less adequate basis 
functions that Coulomb Sturmians, because only the Sturmians possess the correct nuclear cusp 
and radial behavior. 

4 Programming strategy 

Firstly, the ideal ab initio code would rapidly switch from one type of basis function to another. Sec- 
ondly, the chemistry of molecular electronic structure must be used to the very fullest extent. This 
implies using atoms in molecules (AIM) and diatomics in molecules (DIM) at the outset, following 
Bader (in an implementation due to Rico et al [50] and Tully [H] implemented in our previous 
work [55], respectively. The natural choice of atomic orbitals, i.e. the Sturmians or hydrogen-like 
orbitals lend themselves to the AIM approach. To a good approximation, core eigenfunctions for 
the atomic hamiltonian remain unchanged in the molecule. Otherwise, atom pairs are the natural 
choice, particularly if the Coulomb resolution recently advocated by Gill is used. This leads us to 
products of auxiliary overlaps which are either literally one- or two- centered, or have one factor 
of the product where a simple potential function needs to be translated to one atomic center. The 
Slater basis set nightmare of the Gegenbauer addition theorem is completely avoided. Naturally, 
the series of products required for, say a four-center two-electron integral may require 10 or even 20 
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terms to converge to chemical accuracy, when at least one atom pair is bound but the auxiliaries 
are easy to evaluate recursively and re-use. Unbound pairs may be treated using a smaller number 
of terms since the integals can be predicted to be small, using a Schwarz inequality. 
Now, the proposed switch in basis set may also be accomplished just by re-evaluating the auxiliary 
overlaps. Furthermore, the exchange integrals are greatly simplified in that the products of overlaps 
just involve a two-orbital product instead of a homogeneous density. The resulting cpu-time growth 
of the calculation is n 2 for SCF, rather than n 4 . Further gains may be obtained by extending the 
procedure to post-Hartree-Fock techniques involving explicit correlation, since the ri2 _1 integrals 
involving more than two electrons, that previously soon led to bottlenecks, are also just products 
of overlaps. 

5 Atom pairs: solving Poisson's equation 

All the molecular integrals over CS required for standard SCF may be evaluated using analytical 
two-center terms based on the solution of Poisson's equation for the Coulomb potential in an 
ETO basis. This uses the Spectral forms (involving incomplete gamma functions and regular and 
irregular solid harmonics) defined initially in [52[ (63[ [53] and subsequently generalized to ensure 
numerical stability as shown in a brief summary below. 

Recalling the definition of a Slater Type Orbital: 



Then, (from the Spectral forms in [63]), the potential due to this distribution is immediately written: 



Where g is short for g(r) and F(r) is given below, with a suitable variable of integration; u: 



This expression is used to write all radially dependent one and two-center integrals in analytical 
closed form. 

The next section describes a more profound advance, that reduces the atom-pair evaluation to one- 
electron overlap-like integrals. It is related to the Poisson equation technique, as detailed in [63] 
and [33]. 

6 Avoiding ETO translations for two-electron integrals over 3 and 
4 centers 

Previous work on separation of integration variables is difficult to apply, in contrast to the case for 
gaussians [33] cf [ID]. Recent work by Gill et al [15] proposes a resolution of the Coulomb operator, 
in terms of potential functions (pi, which are characterized by examining Poisson's equation. In 
addition, they must ensure rapid convergence of the implied sum in the resulting expression for 
Coulomb integrals J±2 as products of "auxiliaries" i.e. overlap integrals, as detailed in [15] . 



Xn, I, m, c( r >0,0) = r 
Define the radial factor g(r) : 
g {r) = r n - l e-^ r . 



11,(5) = r 2 F(r) , 




l-l 
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This is based on separating the variables of by determining suitable functions of r\ and r 2 that 
treat these variables equivalently and constitute a complete set which orthogonalises the Coulomb 
operator. The associated potentials fa provide an expansion, or resolution of — similar to that of 
the identity (using the summation convention): 

\9i >< 9i\=I (6.3) 

This is the completeness property of a set of orthonormal functions, within a particular Hilbert 
space. Similarly, for the Coulomb operator, suitable potentials give: 

— = | fa >< fa | (6.4) 

This Coulomb resolution is based on a complete set of functions which may be determined such 
that they impose the identity as matrix representation of the Coulomb operator, — in this basis: 

< fi\ — \fj > = <% (6-5) 

T\1 

The completeness relation for the associated potentials can also be written in the form (8.19) The 
functional expression of the above gives: 

— = fa(ri)fa(r 2 ) (6.6) 

The potential functions fa, are solutions of Poisson's equation. The functions chosen may also be 
based on Coulomb Sturmians (see the work by Avery). 

Completeness of the functions fi allows us to expand a density in terms of them: 

< p(r) | =< p(r) \ — \fi(r) > < fi(r)\ (6.7) 
J is re-written: 

J12 = < p(n) | — \p(r 2 ) >= < p(n) | — |/i(n) >< fi(n)\ — \ fj(r 2 ) >< fj(r 2 )\ — \ p(r 2 ) > 
ri 2 ri 2 r 12 r 12 

(6.8) 

summing over i and j. Introducing the orthogonalised operator from 8.20 to resolve the two- electron 
integral into a sum of products of 1-electron overlap-like integrals: 

J12 = < p(^i) I — \fi( r i) > < fi{ r 2)\ — I pi r 2) > with implied summation over i (6.9) 
rn r 12 

And recalling the defining relation for potentials (i.e. one electron functions of a single radial 
variable) : 

\ — \fi(r) >= \fa{r) > (6.10) 
ri2 



J\ 2 = < p{ri) 4>i{r\) > < 4>i{r 2 ) p(r 2 ) > with implied sumation over i (6.11) 
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This technique can be readily generalized to exchange and multi-center two-electron integrals. 
Note, however, that the origin of one of the potential functions only may be chosen to coincide with 
an atomic (nuclear) position. 

Define potential functions (pi in the scope of a Coulomb operator resolution, as follows: 

r+oo 

4>ni{r) = I h n (x)ji(rx)dx with ji(x) denoting the spherical Bessel function (6.12) 
Jo 

Here, h n (x) is the n th member of any set of functions that are complete and orthonormal on 
the interval [0, +oo), such as the n th order polynomial function (i.e. polynomial factor of an 
exponential). The choice made in [15] is to use parabolic cylinder functions (see also another 
application |105j ). i.e. functions with the even order Hermite polynomials as a factor. This is not 
the only possibility and a more natural and convenient choice is based on the Laguerre polynomials 
L n (x): Define: 

h n (x) = V2 L n (2 x)e~ x (6.13) 

These polynomial functions are easy to use and lead to the following analytical expressions for the 
first two terms in the potential defined in (6.2): 

, , /— tan _1 (r) 

V 00 {r) = V2 K -l (6.14) 



r- r tan 1 (r) 2 , 

Vio(r) = V2 [ — li - irT - fy ] (6.15) 

Furthermore, higher n expressions of V n o(r) all resemble (6.5) (see |16] eq (23)): 

Vno{r) = V§ V t(-D fc Sin(2fct r" 1(r)) ) (6-16) 
r i K 

and analytical expressions of V n i(r) with non-zero I are also readily obtained by recurrence. 
The auxiliary overlap integrals < p(r\) <f>i(r\) > and < 4>i{ r 2) p{ r 2) > w iH involve densities obtained 
from atomic orbitals centered on two different atoms in most multi-center two-electron integrals. 
The integrals required in an ETO basis are thus of the type: 

< Mn) Mn) Mn) > (6.17) 

Such integrals appear for two-center exchange integrals and all three- and four center integrals. 
Note that exchange integrals require distinct orbitals tfj a and tp^. In the atomic case, they must 
have different values for at least one of n,l,m or £. In the two-center case, the functions centered 
at a and b may be the same. The product does not correspond to a single-center density: it is 
two-centered. The above equation then illustrates the relationship to the one-electron two-center 
overlap integral, although it clearly includes the extra potential term from the Coulomb operator 
resolution. 

The overlap integrals may be evaluated by separating the variables in prolate spheroidal co- 
ordinates, following Mulliken and Roothaan [12] and using recurrence relations in |18j : 



S 



S( ni ,h,m, n 2 , l 2 ,a, p) = ^1+1/2^+1/2 ^ 2ni )\ {2n 2 W 1/2 s(n x hmn 2 l 2 a$) 

= N(ni,n2,af3)s(ni,li,m,n2,h,0(l3) 

where: a = k\R and /3 = k 2 R- The ki, k 2 are Slater exponents. 
The core overlaps are given by: 

s(ni,h,m,n2,h,a,/3) = f f exp |~^( a + P)^ ~ ~ P) v \ + v ) ni {^ ~ v) n2 T({i, v)d^dv 



1 J-i 



+ ?~f> 

r a ~ r b 



v 



R 



r a and r b are the instantaneous position vectors of the electron from the two centers labeled a 
and b, respectively and separated by a distance R. We also define, suing the normalised spherical 
tensors S: 



The core overlaps then take the form: 



s(m,h,m,n2,h,a,p) = D hh)m Y^Y l x j A i |i(a + /3)| Bj ji(a - 



is a matrix with integer elements uniquely determined from n, I and m. It is obtained as a 
generalised binomial coefficient, in the expansion of (r a — r b ) n (r a + r^)™ 

j 2 ,m is a coefficient that is independent of the principal quantum number. It is obtained upon 
expanding the product of two Legendre functions in this co-ordinate system. Symmetry conditions 
imply that only m\ = m 2 = ra lead to non-zero coefficients. 
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Here, recurrence relations on the auxiliary integrals A and B lead to those for the requisite core 
integrals [151 ITT] . 

This assumes tacitly that the potential obtained from the coulomb operator resolution be centered 
on one of the atoms. Whilst this choice can be made for one pair in a four-center product, it cannot 
for the second. There remains a single translation for this potential in one auxiliary of the two in 
a product representing a four-center integral and none otherwise. The structure of these potential 
functions in (6.5) and (6.6) shows that the translation may be accomplished readily in the prolate 
spheroidal co-ordinates. This point is addressed in detail in a submitted manuscript [36J. 
This method obviates the need to evaluate infinite series that arise from the orbital translations 
efficiently. They have been eliminated in the Coulomb operator resolution approach, since only 
orbitals on two centers remain in the one-electron overlap-like auxiliaries. These can be evaluated 
with no orbital translation, in prolate spheroidal co-ordinates, or by Fourier transformation [36] 116]. 

7 Numerical results compared for efficiency 

Consider the H2 molecule and its dimer/agregates. In an s-orbital basis, all two-center integrals are 
known analytically, because they can be integrated by separating the variables in prolate spheroidal 
co-ordinates. A modest s-orbital basis is therefore chosen, simply for the demonstration on a rapid 
calculation, for which some experimental data could be corroborated. 

The purpose of this section is to compare evaluations using the translation of a Slater type orbital 
basis to a single center (STOP) [55] with the Poisson equation solution using a DIM (Diatomics 
in molecules or atom pair) strategy and finally to show that the overlap auxiliary method is by 
far the fastest approach, for a given accuracy (the choice adopted is just six decimals, for reasons 
explained below). 

H2 molecule with interatomic distance of 1.402d0 atomic units (a.u.) The first table (Table 1) 
assembles the full set of all Coulomb integrals; with one and two-centers evaluated using STOP, 
Poisson and overlap methods. Exponents may be found from the atomic integrals which do not 
include the constant factor (5/8 here). 
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Table 1. 



AOs (zeta) 




ls a 2 


2s a \ 


2s a2 


ls i 1.042999 


1.042999 











lsa2 1.599999 


0.934309 


1.599999 








2s a i 1.615000 


0.980141 


0.870304 


1.615000 





2s a2 1.784059 


0.901113 


0.923064 


1.189241 


1.784059 


ls bl 1.042999 


3.455363 


0.364117 


0.659791 


1.621644 


ls b2 1.599999 


0.433097 


0.332887 


0.635867 


1.541858 


2s bl 1.615000 


0.323691 


0.248050 


0.529300 


1.276630 


2s b2 1.784059 


0.402387 


0.324872 


0.636877 


2.014196 



Table 2. 

Atomic exchange integrals (6 distinct single center values between pairs of different AOs). 



AOs (zeta) 


Label 


[a(l)b(2)a'(2)b'(l)j 


value 


comment 


ls al 1.042999 


1 


1212 


0.720716 




ls a2 1.599999 


2 


1313 


0.585172 




2s al 1.615000 


3 


1414 


0.610192 




2s a2 1.784059 


4 


2323 


0.557878 




ls bl 1.042999 


5 


2424 


0.607927 




ls b2 1.599999 


6 


3434 


0.602141 




2s bl 1.615000 


7 


2121 


0.720716 


= 1212 


2s b2 1.784059 


8 


3232 


0.585172 


= 2323 



Table 3. 

Two-center exchange integrals. All pair permutations possible. Some are identical by symmetry. 



Labels 


value 


comment 


1515 


0.319902 




1516 


0.285009 


= 1525 


1517 


0.325644 


= 1535 


1518 


0.324917 


= 1545 


1527 


0.291743 


= 1536 


1528 


0.293736 


= 1547 


1538 


0.329543 


= 1548 


2525 


0.260034 




2516 


0.254814 




2517 


0.290533 




2518 


0.290149 





The two-center integrals are dominated by an exponential of the interatomic distance and thus all 
ave values close to 0.3. The table is not the full set. All '15' terms, involving ls a i(l) ls b i(2) are 
given, to illustrate symmetry relations. 

Note that this is by no means the best possible basis set for H2, since it is limited to I = functions 
(simply to ensure that even the two-center exchange integral has an analytic closed form). 
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The total energy obtained for the isolated H2 molecule is -1.1284436 Ha as compared to a Hartree- 
Fock limit estimate of -1.1336296 Ha. Nevertheless, the Van der Waals well, observed at 6.4 au 
with a depth of 0.057 kcal/mol is quite reasonably reproduced [52"] . 



Dimer geometry: rectangular and planar. Distance between two hydrogen atoms of neighboring 

molecules: 6 au. Largest two-center integral between molecules: 4.162864 10 -5 . (Note that this 

alone justifies the expression dimer-the geometry corresponds to two almost completely separate 

molecules, however, the method is applicable in any geometry). 

Timings on a Power 6 workstation, for the dimer (all 4-center integrals in msec): 

STOP: 12 POISON: 10 OVERLAP: 2. 

Total dimer energy: -2.256998 Ha. This corresponds to a well-depth of 0.069kcal/mol, which may 
be considered reasonable in view of the basis set. The factor limiting precision in this study is the 
accuracy of input. The values of Slater exponents and geometric parameters are required to at 
least the accuracy demanded of the integrals and the fundamental constants are needed to greater 
precision. 

7.1 The NMR nuclear shielding tensor. 

More complete work is referred to here and the present description is a brief summary |87t [88| I9T] . In 
NMR, the nuclear shielding tensor is a second order perturbation energy correction, for derivatives 
with respect to the nuclear dipole moment and the external field. 

The perturbed Fock matrix element when including the effect of the external field contains both 
one and two electron terms. In this example, we focus on the one electron terms. 
The purpose of the present section is to give a case study of one of the contributing energy integrals 
involving the dipole l/rfj operator. 

In the applied magnetic field, the question of gauge invariance must be resolved. A method of 
circumventing the problem was devised by Ditchfield using the London GIAO [92]. These Gauge 
Including Atomic Orbitals reduce to STO for zero field and contain the required phase factor 
otherwise [871 1881 fTOT] . 

The integrals were evaluated for GTO at zero field and nuclear shielding tensor or chemical shifts 
have been available since Gaussian 72 based on this pioneering work [27] and distributed to aca- 
demics by QCPE. It is nevertheless important that users input the appropriate structure in order 
to obtain accurate chemical shifts corresponding to the species studied and note that for work in 
solution (or in solids) some structural changes may occur. 
Define the nuclear shielding tensor as a second order energy perturbation: 




(7.18) 



with /2at the nuclear dipole moment of nucleus N and Bq the external field. |0) is a closed shell 
ground state Slater determinant, a and /3 stand for cartesian coordinates. 



A coupled Hartree-Fock treatment of the above equation leads to j87j [88] 1100] : 
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where P^ and P^ 1 ' 1 ^ are the density matrix of zero order and first order with respect to the 

external magnetic field, ha is the core Hamiltonian of the first order with respect to nuclear 
dipole moment and is the second order one-electron Hamiltonian with respect to the nuclear 

moment fi a and the external field Br. 

The non-zero orders in (|7.19p involve integrals which are absent from ab initio Hartree-Fock calcu- 
lations. In this work, we focus our attention on integrals involving l/r^y in their operator. These 
integrals appearing in the second order expression for the approximate perturbed Hamiltonian: 



V,<tf- 47r2m P \ Xfl 



r v ■ r N d a/ 3 - r^ a r N> p (R^v A r)p{fjs[ A V) 



+ 



' N ' N 



Xu (7.20) 



The integral which we have chosen to investigate in detail within the Fourier transform approach, 
is the three-center one electron integral: 



r v ■ r N 5 a p - r VtCl r N) p 



Xu) (7.21) 



here rjv is the instantaneous position of the electron with respect to the nuclei N. 

Analytical treatment. 

This is to be found in the Appendix. 

The above algorithm is available in Fortran, within the STOP (Slater Type Orbital Package) set 
of programs, at the coupled perturbed Hartree-Fock level with the ETOs axpanded in Slater type 
orbit als. 

DFT coding proves more accurate for NMR chemical shifts because it accounts for the majority 
of the electron correlation energy. In this case, the ETOs are fitted to large Gaussian expansions, 
following the algorithm in |105j and Gaussian03 is subsequently used. 

Application. 

The 15 N chemical shifts measured for a set of benzothiazoles are evaluated with the above expres- 
sions. These molecules possess a ring nitrogen and have been studied previously in our group |101j . 
The measurements were made in natural abundance. The intensity of signals due to the nitro- 
gen must be amplified by a 2-D NMR technique involving cross-polarization to benefit from the 
intensity of proton resonances coupled to that of the 15 N in the molecule. 

The in vivo NMR benefited from measurements by Bruno Combourieu: these molecules are me- 
tabolized by bacteria and researchers in the group try to follow the pathway by NMR. Since such 
studies are very difficult to do, we tried calculating some chemical shifts accurately from structures 
to assign them (see acknowledgements). 

The Y substituent, generally a hydroxide was found to be in the position indicated (for mechanistic 
reasons, it is the only accessible and stabilised position for ring hydroxylation which has been found 
to take place in vivo after experiments in our group). 

In solution, these molecules undergo a tautomeric equilibrium reaction transferring a proton towards 
this nitrogen as shown in the figure below (also used for nomenclature; P= protonated on resonating 
nitrogen). 

Summary of NDDO-PM3 fitted STO molecular-site calculations on unprotonated tautomers (b). 
It is possible to conclude when the Gaussian03/PBE 6-311+- |-G(2d,p) calculation (c) differs sub- 
stantially from the measured value (a) (ppm/C^iVC^) that the resonating nitrogen is mostly 
protonated. This serves as a guideline for ab initio structures studied for these equilibria. 
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molecule 


substituent 


a 


b 


c 


BT:benzothiazole: 


No X 


-72.5 


-71.8 


-61.4 


OBT: 


x=o 


-238.8 


-238.9 


-133.3 


OHOBT: 


X=0; Y=OH 


-240.4 


-239.9 


-135.3 


ABT: 


X=NH 


-153.1 


-152.1 


-131.6 


OHABT: 


X=NH; Y=OH 


-153.1 


-153.6 


-132.3 


MBT: 


x=s 


-199.6 


-199.9 


-79.9 


OHMBT: 


X=S; Y=OH 


-205.5 


-205.5 


-83.2 


MTB: 


X=N(CH 3 )CONHCH 3 


-124.0 


-125.4 


-141.0 



The above results prompted use of a structure, protonated on the resonating N, (denoted P) to 
obtain the zero-order wave-function, in all cases apart from benzothiazole (BT) and ABT. Below, 
the same cases are treated in the DFT work. 

Note that the basis sets including hydrogen-like orbitals perform better than the STO basis sets 
that in turn improve upon dense-core Gaussian basis sets [6-311+- |-G(2d,p)]. 

Next, examining the generalised basis sets, compare b,c, d, and f with the measured chemical shifts 
and evaluate the difference in ppm. Note that 
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TABLE: DFT Calculations. Differences between calculated and observed N chemical shifts for 
commercial benzothiazoles and some metabolites (in ppm). 

a-Measured values with respect to nitromethane standard in deuterated methanol solvent 
(B. Combourieu in |101| ) error bars of 2 ppm. 
b-Coupled perturbed STO. 

c-Gaussian |89] augmented with hydrogen-like AOs (c.f. Coulomb Sturmians a = 1. 
d-Gaussian [89] , 
e-Generalised ETO a = — 1 
f-Generalised ETO a = -2 

Note, b through f involve solvation models, detailed below. 



molecule 


substituent 


a-c 


a-b 


a-d 


a-e 


a-f 


BT:benzothiazole: 


No X 


1.3 


8.3 


11.1 


1.2 


2.8 


POBT: 


x=o 


4.6 


11.7 


20.0 


3.8 


5.3 


POHOBT: 


X=0; Y=OH 


4.5 


7.4 


14.9 


2.9 


5.2 


ABT: 


X=NH 


1.1 


3.8 


21.5 


0.9 


3.1 


OHABT: 


X=NH; Y=OH 


4.5 


10.1 


20.8 


2.8 


6.1 


PMBT: 


x=s 


3.0 


11.2 


21.2 


2.1 


4.5 


POHMBT: 


X=S; Y=OH 


2.5 


10.1 


18.8 


1.7 


4.8 



Basis sets augmented with hydrogen-like orbitals are within 5 ppm of the experimental values 
(measured within 2 ppm) for the discrete solvated model. This model explicitly includes several 
deuterated methanol molecules to cater for the specific hydrogen bonding interactions. 

a Measured chemical shift for ring nitrogen. 

b STO: DFT PBE 6-311++G(2d,p) calculations with two discrete CD 3 OD molecules on OH, NH, 
and SH (one on N, O, S) for minimal total energy. 

c Gaussian 2003 as (b) with hydrogendike orbital DFT PBE aug-6.311G**(2d, p) calculations, 
d Gaussian 2003 DFT PBE 6-311++G(2d,p) calculations. 

The content of this table is original and based on the previous work of the author [36] i.e. geometries 
are re-optimized from the co-ordinates of |36| . 



Conclusion. 

Another step on the way to ab initio ETO basis nuclear screening tensor calculations has been 
accomplished. 

It is essential to use a basis set which comprises orbitals with the correct nuclear cusp behavior. 
This implies a non-zero value of the function at the origin for spherically symmetric cases and 
satisfying Kato's conditions. Hydrogen-like atomic orbital basis sets therefore perform better than 
Slater type orbitals which are an improvement upon even large Gaussian basis sets. 
The NDDO-PM3 molecular site approach has the advantage of rapidity. Calculations take about 
a minute instead of 50-75 hours on the IBM-44P-270. They cannot be systematically improved, 
however once the site Slater exponents have been fitted. Note that the 2s Slater exponent fluctuates 
wildly in fits, providing further evidence that shielding must be of the form (2-r) for the 2s ETO. 
Fundamental work on orbital translation is also in progress to speed up these calculations within 
the test-bed of the STOP programmes |551 155| 153] . 
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The interplay of these discrete molecule solvent models and accurate in vivo NMR measurements is 
satisfactory, in that the structures postulated give calculated chemical shifts to similar accuracy as 
obtained for experimental values (on the order 2ppm). It should be stressed that energy minimiza- 
tion in this case does evidence directional hydrogen bonds but can lead to several possible solvent 
geometries. Further study, using molecular dynamics techniques would be useful in the modeling 
of solvent shells and is planned in the future. In view of the complex systems studied, this is highly 
satisfactory. 



8 Conclusions 

It is a remarkable gain in simplicity that the Coulomb operator resolution [15] now enables the 
exponential type orbital translations to be completely avoided, although some mathematical struc- 
ture has been emerging in the BCLFs used to translate Slater type orbitals and even more in the 
Shibuya-Wulfman matrix used to translate Coulomb Sturmians. 

This breakthrough in the ETO algorithm strategy represents a well-controlled approximation, anal- 
ogous to resolutions of the identity. The convergence has been shown to be rapid in all cases [16] . 
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10 Appendix: analysis of nuclear dipole integrals for NMR in a 
Slater basis 

The operator, {f v • r^S a ^ — r v ^ a r^^) / (r^J, can be expressed as a combination of terms involving 
cartesian co-ordinates. These terms take the following general form: 

3 10-22 

r N 

where X N, i stands for cartesian coordinates of the electron with respect to the nuclei N. 
Now, it is more convenient to express the cartesian co-ordinates as sums of spherical polar co- 
ordinates with their complex conjugates. These co-ordinates are of the general form: 



x 



2tt 




y = ir 
z = r—Y?{0*<p?) 



] [Yf 1 ^-, w ) + Yt{e F , w )) (10 - 23) 



The product of a STO by a cartesian coordinate can be expressed as a combination of STO, since: 
rY/^,^)^^^ =(-l) M £ (lm\L-M\Xm + M) 

A — A m j n ,2 
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(10.24) 



here we have used the Gaunt coefficients [901 1103] and the Condon and Shortley phase convention 
for spherical harmonics YJ' m (#,^, <pf) [86]. 

Consequently, the integral ()7.21|) is reduced to a sum of terms of the form: 



N 



Xv 



(10.25) 



with just a 1/r 2 dependence and where j = —1,0, 1. 



Using the Fourier transform formalism requires the integral representation of the operator involved 
in ([10.251) . We obtain: 



(P) 



2i Yf(0jf,<pjr) 



N 



P 



This immediately allows us to write the inverse Fourier transform: 



't n i ~Tr N 



<Prlr) 



N 



I 

'2^2 



\p\ 



(10.26) 



(10.27) 



Now, this places us in a position to write the Fourier integral for the present term in the NMR 
nuclear shielding tensor calculation. After expanding the Slater type orbitals in terms of B-functions 
and substituting (j!0.27p in (|10.25p . the present integral becomes: 



i 

2^2 



\p\ 



(10.28) 



X TO 1 (Cl > fO|e-*1^(0«,*'- R2)) ( r) dp 

whereas the three-center nuclear attraction integral is : 

Jp-R N 



1 
2^2 



\p\< 



(10.29) 



The three-center dipolar integral (13) appears in a form closely related to that of the three-center 
nuclear attraction integrals required at the HF-SCF level of electronic structure calculation (and 
also used in electronic DFT work). 

In both above integrals note the presence of the common factor in B-function Fourier transform 
work first studied by the Steinborn group: eq (5.4 and 5.5). See also [106] i.e.: 

- ( 10 - 3 °) 

The analytical treatment developed here has not required any hypothesis on the relative posi- 
tion of nucleus (aligned or not) and any restriction on quantum numbers. Consequently, the 
equation (| 10.28[) is completely general and may be directly evaluated from routines available in a 
quantum calculation software. 

Note that such an integral satisfies all applicability conditions of non-linear transformations for ex- 
trapolation described by A. Sidi [97] • Previous work on three-center nuclear integral evaluation [96] 
has been used to develop an efficient program to compute this dipolar l/r% three-center integral. 
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